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Abstract 

Building on a recently improved understanding of the problem of heat flow in general relativity, 
we develop a hydrodynamical model for coupled finite temperature superfluids. The formalism is 
designed with the dynamics of the outer core of a mature neutron star (where superfluid neutrons 
are coupled to a conglomerate of protons and electrons) in mind, but the main ingredients are 
relevant for a range of analogous problems. The entrainment between material fluid components 
(the condensates) and the entropy (the thermal excitations) plays a central role in the development. 
We compare and contrast the new model to previous results in the literature, and provide estimates 
for the relevant entrainment coefficients that should prove useful in future applications. Finally, 
we consider the sound-wave propagation in the system in two simple limits, demonstrating the 
presence of second sound if the temperature is sub-critical, but absence of this phenomenon above 
the critical temperature for superfluidity. 
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I. INTRODUCTION 



A realistic model of a mature neutron star must account for various extreme aspects 
of physics. In addition to the strong gravity, we need to consider the state of matter at 
supranuclear densities. We also need to ask how the composition of matter changes as 
the density increases. Many different options, involving varying levels of exotica, have 
been suggested. For example, at the present time we do not know if the deep neutron 
star core contains deconfined quarks, perhaps in the form of a colour superconductor 
Meanwhile, the composition of the outer core is better understood Just above the 
crust-core transition density (at around 60% of the nuclear saturation density), the matter 
is predominantly made up of neutrons with a small fraction of protons and electrons (to 
make the conglomerate charge neutral). As the density increases muons appear, eventually 
followed by (this is where uncertainties creep in) hyperons and/or quark deconfmement. As 
a result, a neutron star is a bit like a layer cake where the composition of the different 
layers affects the transport coefficients which in turn determine the effect on cooling and the 
viscous damping of oscillations. Key to modelling this physics in a realistic fashion is the 
fact that several regions in the star are likely to exhibit superfluidity/superconductivity. 

If we focus our attention on neutron star oscillations, then the added degree of freedom 
associated with superfluidity leads to the presence of extra families of oscillation modes (3)- 
f| (analogous to the second sound in laboratory superfluids). This aspect of the problem 
has already been considered at some level of detail. However, most previous models ignore 
finite temperature effects. This may seem natural given that a typical neutron star is 
cold (in the sense that the temperature is orders of magnitude below the neutron Fermi 
temperature). While this statement is certainly true, we also need to acknowledge that the 

n 

various superfluid pairing gaps are density dependent [7|. The upshot of this is that there 
should always be regions in a mature neutron star that are near the superfluid transition, 
see Figure HJ In such a region finite temperature effects will dominate. These transition 
regions may be (highly) localised, but it would nevertheless be wise to consider them in 
detail. They may, in fact, have an unexpectedly large influence on the dynamics of the 
system. Consider for example the analogy with the viscous boundary layer at the crust-core 
interface and the suggested effect on the gravitational- wave driven r-mode instability js-ll]. 
In this case it is the finite extent of the boundary layer that leads to it having a large effect 
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FIG. 1: A schematic illustration of the critical temperatures T c for superfluidity of neutrons (singlet: 
solid line, triplet: dotted line) and protons (singlet: dashed line) for representative pairing gaps, 
cf. [7]. Each critical temperature is shown as function of the total (energy) density p. The neutron 
star crust region at lower densities is indicated in grey. Examples of the critical densities at which 
the superfluid transition takes place for a given temperature are indicated as horizontal (orange) 
lines. The regions where thermal excitations are expected to have a considerable effect on the 
dynamics of the neutron component are shown as vertical (darker) bands for a sample temperature 
of 10 9 K. 



on the dynamics. Thus, without considering the problem in more detail we cannot know if 
the superfluid transition region is dynamically important or not. 

In order to address this problem, we need a model that accounts for finite temperature ef- 
fects in neutron star superffuids. Such models have been discussed, in particular by Gusakov 



and collaborators 
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2Hl4|. but the problem needs further attention. Our aim is to develop a 
phenomenological model that (i) enables us to consider finite temperature superfluid dynam- 
ics, and at the same time is (ii) flexible enough to allow for future extensions involving (for 
example) the star's magnetic field and/or the elastic crust. For obvious reasons the discus- 
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sion will focus on the matter components, especially the extra degrees of freedom that enter 
in a superffuid system. Still, one should keep in mind that the complete model involves the 
"live" spacetime of Einstein's theory. We will not discuss the associated (gravitational-wave) 
aspects here, but it is important to remember that the matter dynamics that we consider is 
coupled to the spacetime geometry. In other words, the spacetime metric g a b is not taken to 
be fixed in the following. In any realistic application, one also has to solve (at least a subset 
of) the Einstein field equations. 

As we intend to develop a covariant model, we distinguish between co- and contravariant 
quantities using the space-time metric to relate the two. Spacetime indices are given as 
italic letters from the beginning of the alphabet, a,b,c, .... In cases where the focus is on 
spatial components the corresponding indices are i,j, k, .... We assume that the spacetime 
signature is +2 and the Einstein summation convention applies. The spacetime indices 
should be distinguished from the constituent indices that are used to identify the different 
fluid components in the problem. These are given as roman, generically x, y, z,... and 
specifically (in the present context) n, p, e and s. Repetition of these component labels in a 
mathematical expression does not imply summation. 



II. THE VARIATIONAL APPROACH 



In contrast to the deep core, the conditions in the outer core of a neutron star are relatively 
well understood. In this region a (triplet pairing) neutron superffuid is thought to coexist 
with a (singlet pairing) type II proton superconductor 15[ and a sea of relativistic electrons 
(we will ignore muons in the following, even though they are quite easy to account for). 
In a hot system we also need to account for the entropy. In general, all four components 
can flow relative to one another. A general model would therefore consider four distinct 
fluid components; neutrons (labelled n in the following), protons (p), electrons (e) and the 
entropy (s). In order to focus on the finite temperature effects, we will ignore two of the 
relative degrees of freedom. The first represents the charge current due to a relative flow 
between protons and electrons (as required to maintain the global magnetic field 16]), and 
the second corresponds to the heat flux fis| . In principle, both of these aspects can be 
included in the model (at the cost of making the involved relations less manageable 16|) 
but the added complexity detracts from the key points of the analysis. Hence, we first of all 
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assume that the protons and electrons are locked electromagnetically. In the final model we 
also add the assumption that the entropy is carried along with the only "normal" component 
in the system, the electrons. This leaves us with two fluid degrees of freedom. This should 
be the minimal model for the system of interest. 

We take as our starting point the convective variational principle for multi-fluid systems 



devised by Carter [21J (see [22j for a review). In this model, the fluid dynamics is obtained 
from a matter Lagrangian A, which is taken to be a function of the various scalars that can 
be constructed from the individual fluxes n", where the constituent index x = n, p,e or s. 
Basically, this means that A is a function of 

= ~« , and n% = -n>£ • (1) 

As will become apparent later, the second of these definitions represents effects due to 
the relative flow between the different components. While this flow is generally expected 
to be small in ma gnit ude, its contribution is nevertheless significant since it encodes the 



entrainment effect 
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-25|. This effect relies on quadratic (and higher order) terms in the 
Lagrangian. Hence, it is important that we do not make immediate use of the low relative 
velocities. This assumption should not be incorporated until the relevant momenta have 
been defined. 

Given the general form for A, the variational procedure determines the momenta that are 
canonically conjugate to the fluxes |22|]. We have 



<9A 



dnt 

x y^x 



where 



(4) 

The latter encodes the entrainment between the different components. We will account for 
two distinct entrainment mechanisms. First of all, we know that the entrainment between 
neutrons and protons arises due to the strong interaction. Hence, we will have A pn ^ 0. We 
also know that the entropy entrainment plays a central role in ensuring causality of the heat 



flux 
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18) . Hence, we assume that *4 SX 7^ (with x = n, p), as well. We will not account 



for entrainment between electrons and baryons. 



It is useful to distinguish the entropy flux s a = n" and the associated thermal momen- 



tum 6 a = /4 
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from the other components. This makes the description slightly less 



"symmetric" but it makes it easier to make contact with thermodynamics. We then have 

ds 



Finally, a variation with respect to the spacetime metric tells us that the stress-energy 
tensor for the system takes the form 

n = m+^nx, (6) 

X 

where the generalized pressure is defined as 

* = A-X>X- (7) 

X 

The stress-energy tensor acts as source in Einstein's equations in the usual way. In the 
case where the fluxes are individually conserved, i.e., when the dynamical timescale is much 
faster than that associated with the relevant reactions, we also have 

V a < = , V a s a = . (8) 

In this (the fully variational) case, each component satisfies an equation of motion of form 

2<V [( X, = , (9) 

and 

2s a V [a 9 b] = . (10) 

These Euler equations combine to ensure the conservation of energy and momentum, 
V a T a6 = 0. 



As in the case of relativistic heat flux 17|, [18| , it is useful to focus on a particular fluid 



frame. Basically, we need to decide which observer is measuring the temperature etcetera. 
In the present context, it seems natural to work in the entropy frame (associated with the 
flow of the "normal" fluid). Hence, we introduce an observer with four velocity u a such that 
s a = su a . In order to quantify the relative flow associated with the superfluid components, 
we define the relative velocities i>° as 

< = 7 x n x K + <), nX = 0, 7x = (l-^ r i/ 2; t£=«X. (11) 



In many cases of interest will be small, and it makes sense to work with an expansion 
with these "drift velocities" as small parameters. Rather than working with the individual 
velocities, it is often convenient to introduce the relative velocities to" = vt — v°;. Then the 
entrainment parameter depends on 

< « (l + ^) , (12) 

and we see that the equation of state, A, must contain information at order w^ y in order 
for the momenta, /i*, to be correct at linear order. Problems that require higher order 
information must be considered on a case-by-case basis. For example, in the Hartle-Thorne 



slow-rotation approximation for a two-fluid system 19j (accurate to order w^ y ) one can show 
that contributions originating from order w^ y terms cancel exactly 

In the general case, the number densities measured in the entropy frame (such quantities 
will be denoted by * in the following) are 

n* = -u a n£ = 7 x n x . (13) 

We then have the momenta (associated with the three different fluxes) 

£ = B n n n a + A pn n p a + A sn s a = (B n n* n + A pn n* p + sA sn ) u a + B n n* n v* + A pn n* p v p , (14) 

= B p n p a + A np nl + A sp s a = (B p n* p + A np n* n + sA sp ) u a + B p n* p v p + A np n* n v* , (15) 

and 

Q a = B s s a + A ns n n a + A ps n p a = (B s s + A ns n* n + A ps n* p ) u a + A ns n* n v n a + A ps n* p v p . (16) 

Note that we have not yet locked the electrons (which carry entropy and therefore move 
along with the entropy component) to the protons (we still have v p ^ 0). Given these 
expressions we see that the "chemical potentials" (in the entropy frame) are 

^ = -u'fi = B n n* n + A pn n* p + sA sn , (17) 

/i p = -«X = B p n* p + A np n* n + sA sp , (18) 

and 

Me = -liVa = B " n e , (19) 
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while the temperature 9 follows from 

6 = -u a 9 a = B s s + A ns n* n + A ps n* p . (20) 

We can use these expressions to rewrite the stress-energy tensor and the pressure. This 
is not necessary but it clarifies the physical meaning of the variational coefficients. We need 

<X = ~n>n + < (B n Kvl + .4 P X«) , (21) 

n a X = -n> p + n* p (B?n;vl + A^vX) , (22) 

n>l = -nel^e , (23) 

and 

s a 6 a = -s9 . (24) 

Given these results we have 

* = A + </i n + n> p + n e /i e + s9 - [l3 n (n> n ) 2 + 2^ n n>;« + B p (n^ p ) 2 ] . (25) 

This is the (extended) Gibbs relation for the system, accounting for the presence of three 
distinct flows. 

Let us now explicitly lock the protons to the electrons, i.e. assume that v£ = 0. Then we 
are left with 

* = A + n> n + n p /i p + n c /i e + s6 - B n (n* n v n f . (26) 

Finally, in most situations of practical interest the deviation from co-motion will be small. 
When that is the case we can neglect terms of order v^, leaving us with 

* = A + n n /i n + n p /i p + n c /i c + s9 . (27) 

In the following we will assume that these various approximations are valid for the neutron 
star problem. Should we be interested in a more general context, it is straightforward to 
reinstate the redshift factors, etcetera. 

Given the above relations, we can work out the corresponding form for the stress-energy 
tensor. To do this, we need 

Kf4 w (* - A) u a u b + n n fi n v^u b + n n {B n n n + A np n p + A ns s) «X 

X 

= - A) u a u b + n n /i n (v a n u b + v\]u a ) . (28) 
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The final expression for the stress-energy tensor is 

T a b = V5 a b + (* - A) u a u b + n D p n (vy b + vlu a ) = -Au a u b + ±\ +n n p n + v£u a ) , 

(29) 

where we have used the standard projection 

1 _ab = g ab + u a u b _ ^ 

t is worth noting that the relative flux enters T\ in the same way as the heat flux would 



17J, |18|. This is, of course, as expected. 

Noting that the energy density measured by an observer moving along with the entropy- 
proton-electron conglomerate is 

p = u a u b T ah = -A , (31) 



we recognize the Gibbs relation (1271) as the definition of the pressure ^ = P, with 9 = T 
the temperature. In other words, in the linear case our final expressions are 

T% = pu a u b + P ±\ +n n p n (v a n u b + vlu a ) , (32) 

and 

P + p = n n p a + n p p p + n e p e + sT . (33) 

The simplified problem has two dynamical degrees of freedom (the remaining fluxes) . To 
determine these we need two momentum equations, e.g. ((91) and ( ITOl) . However, if we want 
to compare to the standard single-fluid problem, then it makes sense to first of all consider 

V b T ab = , (34) 

representing the conservation of total energy and momentum. In addition, we will use the 
Euler equation for the neutrons 

2<V[X] = , (35) 

where, cf., Eq. ( THj) . 

pi = p D u a + B n n n v» . (36) 

So far, the neutron superfluidity is only manifest from the fact that the corresponding 
component is allowed to flow relative to the other fluids in the system, i.e. there is no 
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viscous coupling. This is the generic model one arrives at in situations where quantized 
vortices are present in the system. The macroscopic equations are obtained by averaging 
over the vortices. In order to make contact with the standard analysis of laboratory systems, 
it is useful to consider the fact that a "pure" superffuid should be irrotational. In this case 
the flow is potential, such that we have the momentum 

yU« = m nV a n , (37) 

for some scalar potential <j) Q (proportional to the order parameter of the condensate). In this 
case fl35l) is automatically satisfied. Later, we need to make contact with the "orthodox" 
model for superfluidity, centered on the notion of a superffuid (spatial) "velocity" V^. This 



221 ]. and (in the 



quantity represents the scaled momentum of the superfluid component 
entropy frame) we have 

= — V> n . (38) 
m n 

From (I36I) we then see that we must have 

/i n = -«X = -m n u a V a (j) n , (39) 

and 

v: = ± ±\ ^ = ^< . (40) 

m n m n 

This shows that the "velocity" V® is parallel, but not identical, to the relative velocity f". 
As we will see later, the difference is associated with the entrainment effect (arising when 
the neutron effective mass differs from the bare mass). 

III. NEWTONIAN CORRESPONDENCE AND EFFECTIVE MASSES 

The model we have developed is complete, at least formally. However, in neutron star 
modelling there tends to be a void between theory formulations and applications. The main 
reason for this is that the true supranuclear equation of state remains elusive. To date, most 
discussions have focussed either on bulk properties (as required to construct stationary self- 
gravitating models, leading to mass and radius Q) or the transport coefficients associated 



with cooling 26( (reaction rates etcetera). Comparatively little attention has been given to 



the coefficients required to model superfluid dynamics, like the entrainment or additional 



dissipation coefficients 



27|, 



28]. The notable exceptions are the pioneering work of Alpar 
10 



et al 
;he re 



29j and Borumand et al 30] at the Newtonian level and Comer and Joynt 3l| in 
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ativistic case. The most detailed recent work is that of Chamel anc. 



34| for the neutron star crust and Gusakov et al for hyperon cores 



35 



collaborators 



3q . It is also 



worth noting recent discussions of the additional bulk viscosities that arise in superfluid 



systems 



37 



39j . Another important point concerns the fact that many quantum equation 



of state calculations are carried out in a non-relativistic framework. This is natural since 
;he nuclear physics scale is always small enough that one can work in a local inertial frame 



40f| . However, it means that we need to be able to implement non-relativistic equation of 



state parameters in our fully relativistic framework for superfluid dynamics if we want to 
make progress (given the available information). 

To discuss this issue, it is useful to first focus on a zero temperature model. Then we 
are dealing only with two components and many of the issues become conceptually clearer. 
Moreover, as a first stab at the finite temperature problem it makes sense to build on the 
cold model, simply assuming that the thermal effects (say) the strong interaction are small. 
We develop this strategy through this section and then return to the finite temperature 
problem. 

Let us consider the low-relative velocity version of the variational model. Letting the two 
components be the neutrons and the "protons" (n and p, respectively, with the electrons 
locked to the protons) and writing the variation of A in such a way that the relative velocity 



w 



w 2 is explicit, we have [making use of (JI2 



dA 



w 



n n B n + n p ( 1 + — J A 



l np dn„ 



dn r 



^A^dw 2 . (41) 
2cr 



As we want to consider the Newtonian limit of this expression, we have reinstated the speed 
of light (which is taken to be unity in the rest of the discussion). In considering the limit 
c — > oo we use units such that the time coordinate is x° = ct and the spatial components 
of the four velocity u l x — > v l Jc. With these conventions it is still the case that n 2 . = —n^n* 
has the dimension of a number density. The analysis is also simplified by the fact that we 
are already considering low relative velocities (the detailed calculation has been carried out 
in, for example, [22j and we refer the interested reader to that discussion for more details). 
The key point is that the total energy of the relativistic system includes the rest-mass 
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contributions, but these are not present at the Newtonian level. In effect, we have 



A = - ™xn s c 2 - C , (42) 



where L represents the relevant Newtonian Lagrangian [27j (and the relative sign is conven- 



tion). The Newtonian momenta follow from 27] 

' d£ 
dn 



= KX X + K**n\ , (43) 



x/ n x 



where we have expressed £ in terms of the mobility matrix introduced by Chamel 
is, we define K XJ such that 



45|]. That 



C = ~ Hu^">i + W("x) x = n, p , (44) 

where the fluxes n l x = n x v x are (obviously) spatial (relative to u a ) and U depends on the 
internal energy energy. 

Given these definitions we return to f l4"Tj) . and find that 



Hence, it makes sense to identify 
Similarly, we get 

S) ^ = " n " e ° - np ( J + S - 4 " P = ~ m - c2 ~ "° ' (47) 

where we have defined the Newtonian chemical potential as 

At the level of linear relative velocities, we identify 

n n £ n = (m n -VT n P)c 2 + /^. (49) 
To leading order, the relativistic momentum leads to 

(i n = n n B n + n p A np = m n c 2 + ^ , (50) 
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and 

< = (m n - n p K np ) c< + n p K np cv P = c V \ , (51) 
(the factor of c is required to balance the fact that x° = ct), from which we see that 

n n K nn = m n - n p K np . (52) 

The entrainment has an intuitive interpretation in terms of an effective (dynamical) 
mass. This quantity follows by considering the neutron momentum in the proton frame 
(where v p = 0), and vice versa [25] . That is, we have 

m* = m n — n p K np — > n p K np = m n — m* , (53) 

which means that 

n n K nn = m* n . (54) 
A similar analysis for the protons leads to 

n p K pp = ml , (55) 

and 

n n K np = m p -m* p . (56) 
From this we deduce that that the two effective masses are related: 

m n -m* n = — (m p - m*) . (57) 
n Q 

Before we proceed, it is worth noting that we can (obviously) define an analogous effective 
mass m* in the relativistic model. In this case, the argument leads to [using ( )50l) ] 

m* = n n B n = /i n + (m* - m n ) (58) 

This relation shows that, in the absence of interactions (read, entrainment), the effective 
relativistic mass is equal to the chemical potential. This is as expected. 

Summarising the discussion, we have shown how a relativistic model can be "inferred" 
from the Newtonian limit in the sense that we can obtain B x and A xy from the rest masses 
m x , the effective masses m* and the chemical potentials /i^ (which obviously follow from 
the energy functional IX). The procedure is not "unique", but we have to resort to it in lieu 
of relativistic equation of state models. In particular, the analysis provides the translation 
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between the variational model and the non-relativistic microphysics, e.g. the information 
encoded in the entrainment. The relevant information is expressed in terms of a single 
parameter, which can be taken to be the effective mass ml. 

We now want to build on this strategy, extending it to account for finite temperature ef- 
fects. To do this, we need to make contact with the "orthodox" framework for superfluidity, 
because the problem has so far only been considered in that context [l2, 41]. In other words, 
we want to compare our results to ones based on splitting each condensate into a "normal" 
part and a "superfluid" part (each expressed in terms of a mass density). This language 
provides a convenient description of phenomena observed in laboratory measurements. The 
variational model contains exactly the same information 44(, although the system is de- 
composed in a different fashion. The preference of one model over the other is essentially 
a matter of personal choice. After all, one cannot truly decouple the involved degrees of 
freedom (e.g. the entropy and the nucleons, or the superfluid fraction and the condensate). 
They are inextricably linked. 

As already mentioned, the orthodox model focuses on the superfluid velocities (the 
rescaled momenta). In the Newtonian case, we have = p*/m x which leads to a total 
momentum flux 

f x = p x < = Pxx V: + Pxy V; , y + x , (59) 
where the entrainment is now expressed in the elements of the mass density matrix p xy 



23 



241 ]. Comparing to the variational model, we see that (after a bit of algebra and setting 



m n 



m. 



m for simplicity) 



K 1 



det p 



K pp 



m 2 p mi 
det p 



and 



K np 



m 2 p np 
det p 



where 



det p = pppPn 



P 



np 



(60) 



(61) 



For future reference, it is worth noting that we can relate the mass density matrix to the 
effective masses as well. We then have 

-l 



p n m 



Ppp = Pp 



n p l *\ 
m (m — m ) 

n n 1 



p l ^ 
— (m — m ) 

n p (m — m J 



Pn 



Pp 



m 
m* 



(62) 
(63) 
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and 

Pn P = -p P {m - m* 



m [m — m ) 

p n n v 



-Pp ( ) > (64) 



where the approximations are valid whenever n p /n n <C 1. That is, they should be accurate 
for neutron star cores where the proton fraction tends to be below 10%. We will make use 
of these approximations in our discussion below. 

IV. THE ENTROPY ENTRAINMENT 

So far, most of the literature on superfluid neutron star dynamics has ignored finite 
temperature effects. The exception is a series of paper by Gusakov and collaborators [l2|- 



14j [NEW?] .Their approach builds on the orthodox view of superfluidity and makes use of 



the relativistic hydrodynamics prescription of Son |46] . We want to make contact with this 
model in order to infer the key finite tem per ature effects. In particular, we need to be able 
to use the results of Gusakov & Haensel 4l|, who express the finite temperature effects in 
terms of the mass density matrix. In order to make use of these results we first of all need 
to translate them to the non-relativistic limit of our framework, and then elevate the model 
to relativity. This will not lead to a "complete" model, but it is the best that we can do 
given the lack of microphysics information relating to the problem we are considering. 

The starting point for the discussion is a three-component model, which accounts for the 
superfluidity of both neutrons and protons, and where the thermal excitations are expressed 
in the form of quasi-particles (index qp below). Then the total non-relativistic fluxes can 
be expressed as; 

Jn = Am V£ + Pnp V£ + (p n ~ Pnn ~ Pnp) v\ p , (65) 
fp = Ppp VI + p pn V* + (Pp - P PP - Ppn) < p , (66) 

(here, and in the following, we adapt the notation in such a way that the comparison with 
the variational description in Section II becomes transparent). These expressions account 
for the presence of entrainment through the p xy coefficients in the manner discussed in the 
previous section. It has also been assumed that the "normal velocities" of protons and 
neutrons are identical, represented by v l qp (this is a true transport velocity, in contrast to 
V£). The quasiparticles carry the entropy in the model. 
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The flux expressions are quite intuitive. The density of neutrons flowing with V£ is p nx 
which leaves p n — p nn — p np to flow with v qp (and similarly for the protons). The quasiparticles, 
representing the normal fluid, are only present at finite temperatures. In the T = limit we 
have 

Pn = Pnn + Pnp , (67) 
Pp = Ppp + Pnp , (68) 

bringing us back to the expressions discussed in the previous section. At a finite temperature 
these relations no longer hold, the difference being the presence of thermal excitations. 
This effect was quantified by Gusakov and Haensel 4l|, working in the normal-fluid frame. 
Basically, this means that we let v qp = which leads to an entrainment matrix with the same 
functional form as in the T = case. The similarity is, however, deceptive as the coefficients 
are now temperature dependent. In particular, the three coefficients p xy are independent, 
so we have a three-parameter entrainment model. 

In order to compare the variational model to this description, it is useful to first consider 
a general frame not tied to any of the components. Identifying the quasiparticles with the 
entropy carrying component (the normal fluid) we have v qp = v l a . Assuming that all relative 
velocities are small enough that we can linearise, the momenta will be given by 

pi = n n K nn vi + n p K np v P + sA sn vl , (69) 

and 

pt = n n k np vi + n p K pp v l p + sA sp vl . (70) 



Here we have introduced the entropy entrainment [171. l47j]. encoded in ^4 SX . We have assumed 
that this quantity is trivially related to the corresponding relativistic parameter (this makes 
sense since A np = K np in the cold model). In general, the coefficients of the mobility matrix 
will be affected by temperature, so it will not be the case that K xy = K xy . Of course, this 
depends on the physics involved. One may, for example, assume that the standard entrain- 
ment between neutrons and protons will be essentially the same as in the cold model because 
it originates from the strong interaction, which is only weakly temperature dependent. The 
coupling to the quasiparticles (represented by the entropy carrying component above) com- 
plicates the situation. We need an argument that links the three parameters K np , A sn and 
A sp to the p xy elements of the mass density matrix. 
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If we want to compare (|69|) and (ITUl) to the previous relations, fl65|) and (I66I) . we need to 
invert them to get the respective fluxes. Thus, we arrive at 

777 J*C PP 777 hC np 777 / ~ - \ 

£ = mn n v l n = ^p l n ^p* + =- ( A sp K np - K pp A sn ) svl , (71) 

n detK n detK p det ^ V / s ' v ' 



where 



We also have 



det K = K nn K pp - ^K np y . (72) 



>77 /V^ nn TC n P 777 / „ \ 

jl = mn p < = ^p\ ^ + =r ( ^ sn K np - >Y n M sp ) st£ . (73) 

p Pp detiT P detiT det if V 7 s 

If we introduce the superfluid "velocities" then these relations take the form; 

m 2 K pp m 2 k np 777 / ~ \ 

f n = "-^Jl^v: - —^V* + — Ur ( A sp K np - K pp A sn ) svl , (74) 
det K 11 det K p det K V V s ' 

and 

m 2 K nn m 2 K np ?77 / ~ \ 

= !iiil_v? - ^^V n * + ( A sn K np - K™A sp ) svl . (75) 

p det if p det if det if V J 

At this point, the correspondence with the fluxes in the orthodox formulation is obvious. In 
particular, we see that the quasiparticle contribution vanishes as s — > 0, as expected. It is 
also straightforward to establish that (1601) still holds, although now for the K xy coefficients. 
By comparing to ( 1651) and ( 1661) . and solving for the two entropy entrainment coefficients we 
find 

in 

= sdet p ^ det 9 ~ PppPn + PnpPp> ' ^ 

and 

777 

*4 sp = , , (det p - p nn p p + p n pp n ) . (77) 
sdet p 

This essentially completes the translation between the two models. In particular, we see 
how the entropy entrainment coefficients are related to the mass density matrix elements 
from Gusakov and Haensel j^]. That is, we have all information required for a variational 
finite temperature model for superfluid neutron stars. 

In order to use these results to create a workable relativistic model, it first of all makes 
sense to assume that .4 SX is not affected by relativistic effects (which seems reasonable given 
the discussion in the previous section). Next, we extend the analysis from Section III to 
include the entropy component. This leads to the chemical potential being given by 

p D = n n B n + n p A np + sA sn = m n c 2 + < , (78) 
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and (to leading order in the c — > oo limit) 



n n B n = n n K nn . (79) 

As in the cold case this identification (together with the analogous result for the protons) 
allows us to infer the parameters required in the relativistic case. The model is not unique, 
but it has the correct Newtonian limit by construction. 



V. THE FINITE TEMPERATURE ENTRAINMENT MATRIX 

We have shown how we can infer the entropy entrainment needed to complete the vari- 
ational model from the temperature dependent mass density matrix. This means that we 
can draw on the results of Gusakov and Haensel Q|. They show that the (non-relativistic) 
matrix elements take the form; 

p nn (n,T) = [l-f n (T)]p nn , (80) 
p pp (n,T) = [l-f p (T)]p pp: (81) 

and 

p np (n, T) = [l- / n (T)] [1 - / p (T)] p np , (82) 

^he main temperature behaviour is 



481 ] . These functions vanish in the 



where p xy is only weakly temperature dependent |4i 
encoded in the functions / x , which can be taken from 
zero temperature limit and approach unity as one approaches the critical temperature for 
superfluidity, when T — > T cx . This obviously means that p xn — > as T — > T cn , reflecting the 
fact that all neutrons are "normal" in this limit (similarly for the protons). 

Given this, we may consider the behaviour of the coefficients in the variational model as 
one of the critical temperatures is approached. For this exercise, let us assume that T cn < T cp 
(this would be true near the peak of the triple gap in the example illustrated in Figure []]) . 
Introducing the shorthand notation 

R = PnnPpp - (1 - /n)(l - /p)Pnp > (83) 

we have 

K nn = - 

[1 " fn)R 
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K mi = ( f P F^ , (84) 



m 2 p nn 



and 



K ap = - 



R 



(85) 



(86) 



showing that K xp is regular as T — > TL but f^ nn (and hence B n in the relativistic model) 

nn 

diverges. This behaviour is in sharp contrast with the model of Gusakov et al [121, |4]J , where 
all parameters remain finite. 

However, the divergence of K nn has a natural explanation in terms of the entrainment. 
As the system approaches the critical temperature, the dynamics will become dominated 
by the thermal excitations (an increasing fraction of the neutrons are "normal"). Expressed 
in terms of the entrainment, this has the effect that it is increasingly difficult to move 
the neutrons relative to the entropy (the normal fluid). Essentially, their effective thermal 
mass increases. To see this, let us assume that the entropy is locked to the proton/electron 
component (as in Section II). Then the two momenta take the form 



pi = n n K nn vl + (n p K* p + sA sA j v l p 



(87) 



and 



p l p = n Q K np v l n + \n p K pp + sA sp ) v' p , (88) 
and we see that we can define effective masses m* that include finite temperature effects; 



< = n n K™ , 



(89) 



and 



m* p = n p K pp + s^ sp , (90) 

The first of these clearly diverges as the critical temperature T cn is approached. The second, 
however, remains finite. This becomes clear once we note that K pp is regular and 

(1 - in)PnpPp 



and 



^sn _ ^ 

S 



s 



1 - 



PppPn " 

(1 " fn)R 

PnnPp - (1 - /p)PnpP D 



(i - U)R 

where the former diverges as T — > T cn but the latter remains regular. 



(91) 



(92) 
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As a result of this behaviour, the neutrons lock to the entropy (the "normal" component). 
That is, we must have t>* — > u* in order to ensure that the dynamics remain regular. This is 
apparent from the fact that the neutron momentum (1871) takes the form 



Pn -> 



Pn 



m 



mp p p np 

PnnPpp 



(93) 



(l-A)Pr 

as T — >■ T cn . Demanding that the momentum remains regular simply requires the neutrons to 
lock to the entropy+protons in the limit. This reflects the anticipated physics at the critical 
temperature, and a reduction in the number of distinct dynamical degrees of freedom. 



VI. INSIGHTS INTO THE DYNAMICS 



In order to explore the main features of the finite temperature superfluid model that 
we have developed, let us consider a simplified dynamical problem. The main aim of this 
exercise is to shed light on the transition as the critical temperature is approached. 

We take as our starting point the two-fluid model that results from locking the protons 
(electromagnetically) to the electrons and the entropy, cf., Section II. Meanwhile, the neu- 
trons are superfluid and may move relative to the other components. We consider linear 
deviations from a uniform equilibrium where all fluids move together (the simplest situa- 
tion). We assume that the fluids are in dynamical, chemical and thermal equilibrium in this 
state. This means that all background quantities (including the temperature) are constant 
(in time), and we also have /z n = /i p + /i e . The equilibrium is associated with a frame such 
that u a = [1, 0, 0, 0]. Perturbing this system, the normal fluid will move with u a + 5u a , where 
u a Su a = 0, and the superfluid neutrons move with relative velocity 8v®, where u a Sv^ = 0. 
For simplicity, we will work in the local inertial frame associated with u a . That is, we assume 
that spacetime is flat, which means that all derivatives are partial. 

Let us first consider the conservation law for entropy, which leads to 

d t 6s + sVM = , (94) 

where we have defined the co-moving time-derivative dt = u a \7 a and we have used the fact 
that the perturbed velocities are spatial (i = 1 — 3). Next we perturb the stress-energy 
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tensor and work out its divergence. This leads to 



V a 5T ab = (P + p)u b V a 5u a + (P + p)d t 5u b + u b d t (5P + 8p) 

+ V b 5P + n n p n u b V a 5v^ + n n p n d t 5v b n = . (95) 

Contracting with Ub we arrive at 

d t 8p + (P + p)VM + n^VMn = • (96) 
Meanwhile, we have from the projection orthogonal to u b ; 

(P + p)d t 8ui + Vi5P + n n p n d t 5v n = . (97) 

In these two equations we recognize the standard perfect fluid terms, but also see the influ- 
ence of the superfluid through the terms involving 8v n . Next consider the neutron conser- 
vation law, which leads to 

d t 5n n + n n (V l 5u i + V^<) = , (98) 

while the corresponding momentum equation can be written (the component along u b van- 
ishes identically) 

p n d t 5 Ui + n n B n d t 5vf + V t 5p n = . (99) 

Based on these equations, it is easy to show that only the longitudinal (in the spatial 
sense) degrees of freedom are non-trivial. Hence, we take the divergence of ( 1971) and ( 1991 
in order to obtain scalar equations. Then we can use (1961) and (|98|) to remove the explicit 
dependence on the velocities, ending up with equations relating perturbed scalar variables. 
The relevant equations take the form; 



d 2 t 5p - V 2 5P = , 



(100) 



and 



pI - (p + p)B- 

(P + P) 



d 2 5n n — V 2 5/x n 



p n - n n B n 
(P + P) 



d 2 Sp = . 



(101) 



The first is, obviously, the usual wave equation for sound and the second equation is similar 
(even though the coefficients are a bit messier). If we consider the case where the neutrons 
are not interacting with the other particles (no entrainment), then we have p n = n n B n and 
the second equation collapses to 



d 2 5n n -^(^ )V 2 5n n = . 



Pn \dn n 



(102) 



21 



i.e., a decoupled wave equation with wave-speed [22 ] 



2 n n ( dfi n 



c ' = t \£ ' (103) 



em, leading to the existence of unique "superfluid" oscillations 



The general relations provide a neat demonstration of the existence of two identifiable wave 
degrees of freedom in the prob 
modes in a neutron star core 

Before moving on, it is worth noting that the entropy is slaved to the other perturbations. 
Solving fl96|) and (J9"5j) for the divergence of 5u l and using the result in ({94"]) we see that (after 
time integration) 

5s = — — {Sp - fi n 5n n ) . (104) 

P + p - 72 n /i n 

In order to progress further, we need to introduce a three-parameter equation of state. 
Opting to work with the density, neutron number density and entropy density as our primary 
variables, we have 



and 

In the following we will not state the variables that are held fixed explicitly, for the sake of 
clarity and brevity. In these relations we can use (11041) to replace 5s, effectively reducing the 
problem to one with two parameters. Having done this, we can use the obtained expressions 
in (|97|) and (199]) . leading to two coupled wave equations for 5p and 5n n . 
These equations are rather messy, with the first taking the form 

d 2 t 5p - aV 2 5p - bV 2 5n n = , (107) 

where 

«=( , . , 

dp J P + p — n n p, Q \ds J 



0P ) + ^ (T) • (108) 



and 



i={ ^) do9) 

\dn n J P + p- n n p n \ ds J 



The second equation can be written 

cd 2 5n n - dV 2 5n n - eV 2 5p = , (110) 
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with 



pI-{p + pW 

n n p n -(P + p) 



dp n 



Sp n 



P + p - n n p D V ds 



+ 



and 



p n - n n B n 



dp n 



+ 



p n ~ n n B n 
n n Pn ~ (P + p). 

sp n ( dp n 



'111) 



;ii2) 



;ii3) 



n D p n — (P + p) J \ dp J P + p- n n p n \ ds 
Given these expressions it is straightforward to work out the dispersion relation for plane 
waves in the system. Assuming that the perturbations behave as exp(iut + ikjx^) and 
defining the phase- velocity a = u/k we find that 



(a 2 - a) [co 2 - d) - be = . 



(114) 



It is obviously easy to obtain the characteristic wave speeds in the system, once we introduce 
an explicit equation of state. We will, however, not do this. Instead, let us focus on two 
limiting cases. Consider first the low-temperature limit, i.e. when the system far below the 
critical temperature, deep into the superfluid region. Then we can take the s — > limit to 
obtain the dispersion relation; 



a 



dP\ 
~dp~) 



a 



dp n 



n. 



dP 



dp Q 



. 



;ii5) 



/in \dn n J J p n \dn n J \ dp 

As a sanity check of this result, it is worth noting that in the zero temperature limit we can 
rewrite the expression using n n and n p as the independent variables. Then defining 

dp x 



x p x \dn x J ny 

we arrive at the (more symmetric) dispersion relation 

'9p n 



116) 



( 



a 



J [ a ~ C p) " — ( 



Up f dpp 
P P \dn n 



. 



(117) 



This result is identical to that derived in 



22]. 



The other useful limit to consider corresponds to letting the temperature approach the 
critical temperature for the neutrons. We know from the discussion in the previous section 
that the parameters in the model are singular in this limit. In particular, we know that 



B n ~K nn ^oc 



(118) 
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while the chemical potential p n remains finite. In order for the problem to remain regular, all 
terms proportional to B n in the equations must balance. Thus, we find that (jllOp becomes 

- I -^dHn n + bV 2 5n n + aV 2 6p = . (119) 

Adding this to (11071) (and time integrating) we see that 

Sn n = yf~ p Sp ■ (120) 

That is, one of the dynamical degrees of freedom has been quenched. The remaining wave 
equation can be written 



„ 2A \fdP\ s /8P\ n„ fdP 
which leads to the sound speed 



V 2 Sp = , (121) 



To see that this is, indeed, the expected result we need to take a few more steps. First 
we note that the standard single fluid analysis is based on the total baryon number density 
— n n + n p- Secondly, it is customary to work in terms of the entropy per baryon rather 
than the entropy number density. Thus we define s = s/n^. We can then rewrite the 
thermodynamical derivatives in terms of p and s. This leads to the final result 



dp 

which is the standard expression for the (adiabatic) sound speed in the single-fluid case. 

This model problem demonstrates that the finite temperature model we have developed 
connects the two-fluid region with the normal fluid region in a natural fashion. 

VII. BRIEF SUMMARY AND OUTLOOK 

Working in the context of the variational framework for relativistic multi-fluid systems 
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22], we have developed a model that accounts for finite temperature effects in coupled 
superfluid systems. The model builds on recent progress on the analogous problem of heat 
flow 
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18j and is based on treating the thermal excitations (e.g. phonons) in the systems 



as an additional entropy-carrying fluid. The approach breaks down in the zero-temperature 
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limit, where the excitations propagate ballistically rather than collectively, but should be 
relevant for various systems of astrophysical interest. Our model is designed with the dy- 
namics of the outer core of a mature neutron star (where superfluid neutrons are coupled to 
a conglomerate of protons and electrons) in mind, but the main ingredients are relevant for 
a range of analogous problems. 

The developed formalism allows for detailed studies of the dynamics of neutron star super- 
fluids near the critical density at which superfluidity sets in for a given (local) temperature. 
Setting the scene for more detailed work on this problem, we explored the key role of the 
entrainment between the entropy and the superfluid condensates, arguing that an increas- 
ing "thermal effective mass" couples the components as the critical density/temperature is 
approached. We demonstrated this coupling in the specific case of the sound waves in the 
system. This exercise brings out the expected behaviour; Far below the critical temperature 
the systems exhibits second sound, but this phenomenon is quenched as one approaches the 
edge of the superfluid region. This result has implications for future work on the dynamics 
of neutron star superfluids. In particular it improves our understanding of the normal to 
superfluid transition, a key issue for studies of, for example, realistic neutron star seismology. 
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